Fecha de entrega: 4 de Julio a las 23:59 hrs.
Profesor: Jorge G. Amaya A.
Auxiliar: Aldo Gutiérrez Concha.
Ayudantes: Carolina Chiu y Mariano Vazquez.
Autor: Felipe Urrutia Vargas
Repositorio: github.com/furrutiav/nonlinear_optimization_zoutendijk
En este trabajo nos enfocaremos en resolver el problema de optimización no lineal: $$(P) \quad \begin{cases} \min &f (x) \\ s.a. &Ax \leq b \\ &Ex = e \end{cases}$$ en su generalidad.
Para esto, se realizara una implementacion del metodo Método de direcciones admisibles o metodo de Zoutendijk. Algoritmo que sigue los siguientes pasos:
(0) Sean $\epsilon>0$, $k = 0$, $x_0 \in \mathbb{R}^n$ tal que $Ax_0 \leq b$, $Ex_0=e$
(1) Sea la descomposición $A = [A_1,A_2]^T, b =(b_1, b_2)^T$ tal que $A_1x_k = b_1, A_2x_k < b_2$.
(2) Resolver el problema lineal
$$(D_k) \quad \begin{cases} \min &\nabla f (x_k)^T d\\ s.a. &A_1d \leq 0 \\ &Ed = 0\\ &−1 \leq d_j \leq 1, \quad j = 1, . . . , n \end{cases}$$y sea $d_k$ solución de $(D_k)$.
Si $|| \nabla f (x_k)^T d_k || < \epsilon$, entonces parar.
Si no, ir a (3).
(3) Determinar el paso, resolviendo aproximadamente el problema de minimización unidimensional
$$(L) \quad \begin{cases} \min &f (x_k + \lambda d_k)\\ s.a. &\lambda \in [0, \bar{\lambda_k}] \end{cases}$$mediante el método de Armijo.
Se usa
$$ \bar{\lambda_k} = \min \left\lbrace \frac{(b_2 − A_2 x_k)_i}{(A_2 d_k)_i} : (A_2 d_k)_i > 0 \right\rbrace$$y se considera $\bar{\lambda_k} = +\infty$ cuando $(A_2 d_k)_i \leq 0, \forall i$.
Sea $\lambda_k$ la solución del subproblema $(L)$. Hacer:
$$x_k+1 = x_k + \lambda_k d_k,$$$$k ← k + 1,$$e ir a $(1)$.
Esta seccion se divide en los siguientes capitulos:
Los capitulos Paso i corresponden a implementar el i-esimo paso del Pseudo-algoritmo.
Mientras que el capitulo Clase Zoutendijk corresponde al desarrollo de una clase (paradigma de programacion de python) que permite ingresar y resolver cualquier problema de tipo $(P)$.
La siguiente celda contiene todas las librerias utilizadas en este trabajo.
from scipy.optimize import linprog
import numpy as np
!pip install numdifftools
import numdifftools as nd
import re
import pandas as pd
import time
import plotly.express as px
Requirement already satisfied: numdifftools in c:\users\felip\anaconda3\lib\site-packages (0.9.40) Requirement already satisfied: statsmodels>=0.6 in c:\users\felip\anaconda3\lib\site-packages (from numdifftools) (0.13.0) Requirement already satisfied: algopy>=0.4 in c:\users\felip\anaconda3\lib\site-packages (from numdifftools) (0.5.7) Requirement already satisfied: numpy>=1.9 in c:\users\felip\anaconda3\lib\site-packages (from numdifftools) (1.21.5) Requirement already satisfied: scipy>=0.8 in c:\users\felip\anaconda3\lib\site-packages (from numdifftools) (1.7.3) Requirement already satisfied: pandas>=0.25 in c:\users\felip\anaconda3\lib\site-packages (from statsmodels>=0.6->numdifftools) (1.4.1) Requirement already satisfied: patsy>=0.5.2 in c:\users\felip\anaconda3\lib\site-packages (from statsmodels>=0.6->numdifftools) (0.5.2) Requirement already satisfied: pytz>=2020.1 in c:\users\felip\anaconda3\lib\site-packages (from pandas>=0.25->statsmodels>=0.6->numdifftools) (2021.3) Requirement already satisfied: python-dateutil>=2.8.1 in c:\users\felip\anaconda3\lib\site-packages (from pandas>=0.25->statsmodels>=0.6->numdifftools) (2.8.2) Requirement already satisfied: six in c:\users\felip\anaconda3\lib\site-packages (from patsy>=0.5.2->statsmodels>=0.6->numdifftools) (1.16.0)
WARNING: Error parsing requirements for torch: [Errno 2] No such file or directory: 'c:\\users\\felip\\anaconda3\\lib\\site-packages\\torch-1.11.0.dist-info\\METADATA'
Los parametros del paso 0 son $\epsilon>0$, $k = 0$, $x_0 \in \mathbb{R}^n$ tal que $Ax_0 \leq b$, $Ex_0=e$. Notar que $\epsilon$ y $k$ son arbitrarios. Por otro lado, el punto inicial $x_0$ debe de satisfacer las restricciones lineales $Ax\leq b$ y $Ex=e$. Para verificar la factibilidad del punto inicial, se utiliza la siguiente función:
def step0(xk, A, b, E=None, e=None):
"""
Función que evalua la factibilidad del punto xk
en el poliedro definido por {x : Ax<=b, Ex=e}.
Parametros
----------
xk : np.array (n, 1)
Punto a evaluar
A : np.array (m_A, n)
Coef. de las restr. de desigualdad
b : np.array (m_A, 1)
Lado derecho de las restr. de desigualdad
E : default=None; np.array (m_E, n)
Coef. de las restr. de igualdad
e : default=None; np.array (m_E, 1)
Lado derecho de las restr. de igualdad
"""
feasible = np.all(A@xk <= b)
if str(E) != "None":
feasible = feasible and np.all(np.isclose(E@xk, e))
return feasible
Dado un punto factible $x_k$, el objetivo del paso 1 es descomponer las restricciones de desigualdad por componentes activas e inactivas. Es decir, una descomposición de $A = [A_1,A_2]^T$ y $b =(b_1, b_2)^T$ tal que $A_1x_k = b_1, A_2x_k < b_2$. Para resolver esto, se utiliza la siguiente función:
def step1(xk, A, b, verbose=False, k=-1):
"""
Función que descompone las restricciones de
desigualdad (Ax<=b) en aquellas que son activas
(A1 x = b1) y las que NO (A2 x < b2), donde
A = [A1 | A2] y b = [b1 | b2]
Parametros
----------
xk : np.array (n, 1)
Punto a evaluar
A : np.array (m_A, n)
Coef. de las restr. de desigualdad
b : np.array (m_A, 1)
Lado derecho de las restr. de desigualdad
"""
index_A1 = np.isclose(A@xk, b)
index_A2 = ~index_A1
A1 = A[index_A1]
A2 = A[index_A2]
b1 = b[index_A1]
b2 = b[index_A2]
if verbose:
print(f"(step=1, k={k}) A1.index: {index_A1}")
print(f"(step=1, k={k}) A2.index: {index_A2}")
return A1, A2, b1, b2
Dado un punto factible $x_k$ tal que $A_1 x_k = b_1$, el objetivo del paso 2 es resolver el siguiente problema lineal:
$$(D_k) \quad \begin{cases} \min &\nabla f (x_k)^T d\\ s.a. &A_1d \leq 0 \\ &Ed = 0\\ &−1 \leq d_j \leq 1, \quad j = 1, . . . , n \end{cases}$$Para determinar el gradiente de la funcion en el punto $x_k$, utilizamos la libreria $\texttt{numdifftools}$ con la clase $\texttt{Gradient}$. Para leer la documentacion, revisar aqui.
Para resolver $(D_k)$ en python, consideramos la libreria $\texttt{scipy}$ con el modulo $\texttt{optimize}$. Entra las funciones de optimizacion disponibles, utilizaremos la funcion $\texttt{linprog}$ que se encarga de resolver este tipo de problema de programacion lineal. Para mas detalle, revisar la documentacion disponible aqui. Por defecto, utilizaremos el método simplex para resolver cada problema $(D_k)$.
Llamaremos por $d_k$ a la solución de $(D_k)$. Para la siguiente proposicion $| \nabla f (x_k)^T d_k | < \epsilon$ (criterio de parada), avisaremos con una variable $\texttt{stop}$ su valor de verdad (verdadero o falso).
En la celda que sigue, se encuentra una funcion que resuelve $(D_k)$ y entrega tanto la direccion $d_k$ encontrada como el valor de verdad, del criterio de parada, en una variable $\texttt{stop}$:
def step2(obj, A1, E, eps, method="simplex", verbose=False, k=-1):
"""
Función que resuelve el problema lineal (Dk)
definido por:
(Dk) | min <obj, d>
| s.a A1d <= 0
| Ed = 0
con dk la solución optima.
Parametros
----------
obj : np.array (n, 1)
Costos de la funcion lineal objetivo (grad f(xk))
A1 : np.array (m_A1, n)
Coef. de las restr. de desigualdad activas
E : np.array (m_E, n)
Coef. de las restr. de igualdad
eps : float
Tolerancia para chequear que |obj * dk| sea pequeño
"""
lhs_ineq = A1
rhs_ineq = np.zeros(A1.shape[0]) if str(A1) != "None" else None
lhs_eq = E
rhs_eq = np.zeros(E.shape[0]) if str(E) != "None" else None
bnd = [(-1, 1) for _ in range(obj.shape[0])]
opt = linprog(c=obj, A_ub=lhs_ineq, b_ub=rhs_ineq,
A_eq=lhs_eq, b_eq=rhs_eq, bounds=bnd,
method=method)
dk = opt.x
stop = abs(np.dot(obj, dk))<eps
if verbose:
print(f"(step=2, k={k}) dk: {dk}")
return dk, stop
Dada un direccion de descenso $d_k$ en un punto factible $x_k$, el objetivo del paso 3 es determinar el paso $\lambda_k$. Para esto, se propone resolver aproximadamente el siguiente problema de minimización unidimensional:
$$(L) \quad \begin{cases} \min &f (x_k + \lambda d_k)\\ s.a. &\lambda \in [0, \bar{\lambda_k}] \end{cases}$$Para determinar $\bar{\lambda_k}$, se usa
$$ \bar{\lambda_k} = \min \left\lbrace \frac{(b_2 − A_2 x_k)_i}{(A_2 d_k)_i} : (A_2 d_k)_i > 0 \right\rbrace$$y, en caso de que el conjunto sea vacio, se considera $\bar{\lambda_k} = +\infty$. En la implementacion, $+\infty = 1$.
Parte de este trabajo es resolver $(L)$ mediante el método de Armijo. Este metodo consiste en lo siguiente: Se define $\varphi(t) = f (x_k + t d_k)$ y se elige $\sigma \in (0,1)$ de tal forma que la recta de pendiente $\sigma \varphi'(0)$ se utiliza para acotar superiormente a la funcion en una vecindad del cero; y sobre este dominio encontrar el paso mas largo posible. En efecto, encontrar
$$t^* = \max \{t : \varphi(t) \leq \varphi(0) + \sigma \varphi'(0)t\}$$Para resolver este otro problema, consideramos la siguiente aproximacion:
$$t^* = \max \{t : \exists p\in \mathbb{N}, t=ph, \varphi(t) \leq \varphi(0) + \sigma \varphi'(0)t\}$$donde $h>0$ es dado por el usuario. En este trabajo, se definira $h = \frac{\bar{\lambda_k}}{N}$ con $N>0$ grande y dado por el usuario.
Para determinar la derivada de la funcion $\varphi$ en el cero, utilizamos la libreria $\texttt{numdifftools}$ con la clase $\texttt{Gradient}$. Para leer la documentacion, revisar aqui.
Sea $\lambda_k$ la solución del subproblema $(L)$. Con dicha solucion, hay que generar la siguiente iteracion dado el paso y la direccion:
$$x_{k+1} = x_k + \lambda_k d_k,$$Finalmente, en la celda que sigue se encuentra la funcion que resuelve el problema $(L)$, y retorna el nuevo punto $(x_{k+1})$. Adicionalmente, entrega al valor de verdad de la siguiente proposicion $\lambda_k< \epsilon$. Dicha proposicion, es otro criterio de parada que consideramos produente para el algoritmo.
def step3(f, xk, dk, A2, b2, eps, N=10**2, sgm=0.1, verbose=False, k=-1):
"""
Función que resuelve el problema undimensional (L)
definido por:
(L) | min f(xk + lmbd dk)
| s.a lmbd <= lmbdk
| lmbd >= 0
con lmbdk la solución optima.
Parametros
----------
f : function
Funcion objetivo
xk : np.array (n, 1)
Punto actual
dk : np.array (n, 1)
Direccion de descenso
A2 : np.array (m_A2, n)
Coef. de las restr. de desigualdad no-activas
b2 : np.array (m_A2, 1)
Lado derecho de las restr. de desigualdad activas
eps : float
Tolerancia para chequear que lmbdk sea pequeño
N : int
Tamaño de la discretización equi-espaciada del intervalor
[0, bar_lmbdk]
sgm : float en (0, 1)
Pendiente del paso de armijo
"""
# Lambda maximo
dom = A2@dk>0
bar_lmbdk = np.min((b2-A2@xk)[dom]/(A2@dk)[dom]) if np.sum(dom) != 0 else float("inf")
# Paso de armijo
## Discretización del intervalo [0, bar_lmbdk]
t = np.linspace(0, np.min([bar_lmbdk, 1]), N)
h = 1/(N-1)
## Gradiente de phi en cero
d_phi_0 = nd.Gradient(lambda t: f( xk + dk*t ))([0])
## lmbdk optimo
dom = np.where(f( (xk + dk * t[:, None]).T ) <= f(xk) + sgm*d_phi_0*t)[0]
ix_max = np.max(dom)
lmbdk = ix_max * h * bar_lmbdk
# Nuevo paso
xk1 = xk + lmbdk * dk
stop = lmbdk<eps
if verbose:
print(f"(step=3, k={k}) bar_lmbdk: {bar_lmbdk}")
print(f"(step=3, k={k}) lmbdk: {lmbdk}")
return xk1, stop
En esta seccion se realiza un ensamble de los pasos ya definidos. Adicionalmente, se diseña un metodo denominado $\texttt{solve}$ que ejecuta este ensamble dado ciertos parametros iniciales. Dicho metodo, es parte de la clase principal de esta seccion que denominaremos como $\texttt{Zoutendijk}$.
La idea es que esta clase, logre resolver en su generalidad cualquier problema de tipo $(P)$. Para esto, se añade un lector con expresiones regulares para leer las restricciones. De tal forma que dado un string logre generar las metrices de desigualdades e igualdades automaticamente. Asi mismo, la funciones objetivo es escrita en formato $\texttt{python}$. Por ejemplo, la funcion $x^2$ se introduce por $\texttt{x**2}$.
Ahora bien, consideramos los siguientes criterios de parada para el metodo que resuelve los problemas de tipo $(P)$. Estos criterios son:
Cada uno logra evitar posibles situaciones donde realizar otra iteraciones es marginal al resultado final. Se añade el numero maximo de iteraciones para evitar casos patologicos como por ejemplo entrar en un loop que generen iteraciones ciclicas. Por otro lado, detenerse si la siguiente iteracion es infactible, caso que puede ocurrir por errores numericos. La restriccion de $||\nabla f(x_k)^T d_k|| < \epsilon$ es la que viene por defecto con el metodo de $\texttt{Zoutendijk}$. Sin embargo, se considero tambien $\lambda_k < \epsilon$, que se traduce en un paso pequeño.
En lo que sigue se encuentra la clase ya descrita. Por simplicidad, pasaremos directamente a los resultados:
class Zoutendijk(object):
def __init__(self, fun: str, const: list, info: bool =True):
self.fun = fun
self.const = const
self.info = info
self.set_vars()
self.set_f()
self.set_A_b_E_e()
if self.info:
print(f"""-------------- Model Size --------------
Variables : {self.n}
Constraints : {self.m}
Eq. consts. : {self.m_E}
Ineq. consts.: {self.m_A}
------------------------------------------""")
def set_vars(self):
self.vars = set(re.findall("x\d", self.fun))
for c in self.const:
self.vars = self.vars.union(set(re.findall("x\d", c)))
self.vars = sorted(list(self.vars))
self.dic_vars = {var: f"x[{k}]" for k, var in enumerate(self.vars)}
def set_f(self):
self.rename_fun = self.fun
for k, v in self.dic_vars.items():
self.rename_fun = self.rename_fun.replace(k, v)
def set_A_b_E_e(self):
self.n = len(self.dic_vars)
self.m = len(self.const)
self.m_A = sum(int("<" in c or ">" in c) for c in self.const)
self.m_E = self.m-self.m_A
self.A = None
self.E = None
self.b = None
self.e = None
if self.m_A:
self.A = np.zeros((self.m_A, self.n))
self.b = np.zeros(self.m_A)
if self.m_E:
self.E = np.zeros((self.m_E, self.n))
self.e = np.zeros(self.m_E)
i_A = 0
i_E = 0
for c in self.const:
if "<" in c or ">" in c:
sgn = 1 if "<" in c else -1
left, right = c.split("<=") if sgn>0 else c.split(">=")
self.b[i_A] = sgn*float(right.strip())
t_2 = ""
t_1 = ""
t = ""
for t1 in left.strip().split():
if t1 in self.dic_vars.keys():
v = self.dic_vars[t1]
j = int(re.findall("\d", v)[0])
if t_2 == t_1 == t == "":
self.A[i_A, j] = sgn
elif t_2 == t_1 == "" and t == "-":
self.A[i_A, j] = sgn*(-1)
elif t_2 in ["", "+"] and t == "*":
self.A[i_A, j] = sgn*float(t_1)
elif t_2 == "-" and t == "*":
self.A[i_A, j] = sgn*float(t_1)*(-1)
elif t_1 in self.dic_vars.keys() and t == "+":
self.A[i_A, j] = sgn
elif t_1 in self.dic_vars.keys() and t == "-":
self.A[i_A, j] = sgn*(-1)
t_2 = t_1
t_1 = t
t = t1
i_A+=1
else:
left, right = c.split("=")
self.e[i_E] = float(right.strip())
t_2 = ""
t_1 = ""
t = ""
for t1 in left.strip().split():
if t1 in self.dic_vars.keys():
v = self.dic_vars[t1]
j = int(re.findall("\d", v)[0])
if t_2 == t_1 == t == "":
self.E[i_E, j] = 1
elif t_2 == t_1 == "" and t == "-":
self.E[i_E, j] = -1
elif t_2 in ["", "+"] and t == "*":
self.E[i_E, j] = float(t_1)
elif t_2 == "-" and t == "*":
self.E[i_E, j] = float(t_1)*(-1)
elif t_1 in self.dic_vars.keys() and t == "+":
self.E[i_E, j] = 1
elif t_1 in self.dic_vars.keys() and t == "-":
self.E[i_E, j] = -1
t_2 = t_1
t_1 = t
t = t1
i_E+=1
def f(self, x):
return eval(self.rename_fun)
def solve(self, x0, eps=1e-6, N=10**2, sgm=0.5, max_iter=1000, verbose=True):
if self.info:
print(f"""----------------- Params -----------------
Start point : {x0}
epsilon : {eps}
(Armijo) N : {N}
(Armijo) sigma : {sgm}
(Armijo) max. iter.: {max_iter}
------------------------------------------""")
start=time.time()
# Step 0: Inicializar
k = 0
xk = x0
grad = nd.Gradient(self.f)
# Step 0: Check factibilidad
feasible = step0(xk, self.A, self.b, self.E, self.e)
stop = not feasible
# Step 0: Imprimir estado
f_xk = self.f(xk)
if verbose:
print(f"(step=0, k={k}) factible: {feasible}")
print(f"(step=0, k={k}) xk: {xk}")
print(f"(step=0, k={k}) f(xk): {f_xk}\n")
self.list_xk = [xk]
self.list_f_xk = [f_xk]
self.list_feasible = [feasible]
while (not stop) and (k<max_iter):
# Step 1: Descomposicion [A1, A2]
A1, A2, b1, b2 = step1(xk, self.A, self.b, verbose=verbose, k=k)
# Step 2: Resolver problema lineal (Dk)
obj = grad(xk)
dk, stop_grad = step2(obj, A1, self.E, eps, verbose=verbose, k=k)
stop = stop_grad
if not stop:
# Step 3: Resolver problema (L)
xk1, stop_lambda = step3(self.f, xk, dk, A2, b2, eps=eps, N=N, sgm=sgm, verbose=verbose, k=k)
xk = xk1
feasible = step0(xk, self.A, self.b, self.E, self.e)
stop = (not feasible) or stop_lambda
# Step 3: Imprimir estado
f_xk = self.f(xk)
if verbose:
print(f"(step=3, k={k}) factible: {feasible}")
print(f"(step=3, k={k}) xk: {xk}")
print(f"(step=3, k={k}) f(xk): {f_xk}\n")
k+=1
self.list_xk.append(xk)
self.list_f_xk.append(f_xk)
self.list_feasible.append(feasible)
self.dt = time.time()-start
self.stop_criterion = "|f(xk) * dk| < eps" if stop_grad else "lmbdk < eps" if stop_lambda else "Unfeasible" if not feasible else "max. iter." if k>=max_iter else "?"
self.n_iter = k
if self.info:
if k>0:
print(f"""---------------------------------------------------
Solver : Zoutendijk (simplex/Armijo)
Solution time : {self.dt:.3f} ms
Objective : {self.list_f_xk[-1]: .4f}
Stop criterion : {self.stop_criterion}
{"Max. iterations" if self.n_iter == max_iter else "?" if not feasible else "Successful solution"}
---------------------------------------------------""")
if k==0:
print(f"""---------------------------------------------------
Solver : Zoutendijk (simplex/Armijo)
Objective : {self.list_f_xk[-1]: .4f}
Start point : {"Unfeasible" if not feasible else "?"}
---------------------------------------------------""")
def report(self, feasible=False):
df = pd.DataFrame(self.list_xk, columns=self.dic_vars.keys(), index=pd.Series(range(len(self.list_xk)), name="Iter. k"))
df["f(xk)"] = self.list_f_xk
if feasible:
df["Factible"] = self.list_feasible
return df
def opt(self):
return {
"xf": self.list_xk[-1],
"f(xf)": self.list_f_xk[-1]
}
def error(self, sol, val):
if self.info:
print(f"""---------------------------------------------------
Last point : {self.list_xk[-1]}
Solution : {sol}
Objective : {self.list_f_xk[-1]: .6}
Value : {val: .6f}
---------------------------------------------------""")
error_sol = np.linalg.norm(self.list_xk[-1] - sol)
error_val = abs(self.list_f_xk[-1] - val)
return {
"error_sol": error_sol,
"error_val": error_val
}
def status(self):
return {
"dt": self.dt,
"stop_criterion": self.stop_criterion,
"n_iter": self.n_iter,
"feasible": self.list_feasible[-1]
}
En esta seccion se realizara la resoluciones de dos problemas con el metodo de Zoutendijk.
Para comparar nuestros resultados con el optimo real del problema $(P_1)$ consideramos el motor https://www.wolframalpha.com. La instancia de este problema se puede encontrar aqui. Que entrega la siguiente solución:
| x1 | x2 | f(x1, x2) |
|---|---|---|
| 3.858319799174562844633020514 | 0.2125203012381557330504692291 | 46.90291234143656803690204873 |
sol_P1 = np.array([3.858319799174562844633020514, 0.2125203012381557330504692291])
val_P1 = 46.90291234143656803690204873
%%time
problem_1 = Zoutendijk(
fun="8 * ( x1 - 6 ) ** 2 + ( x2 - 2 ) ** 4",
const=[
"- x1 + 2 * x2 <= 4",
"3 * x1 + 2 * x2 <= 12",
"x1 >= 0",
"x2 >= 0"
],
info=True
)
-------------- Model Size -------------- Variables : 2 Constraints : 4 Eq. consts. : 0 Ineq. consts.: 4 ------------------------------------------ Wall time: 0 ns
Este es un problema de dos variables con cuatro restricciones. Todas ellas de desigualdad.
Inicializamos con un punto incial, epsilon pequeño, sigma igual a $0.2$ y un $N=100$. Donde $N$ es tal que $h_k = \frac{\bar{\lambda}_k}{N}$ es el $h$ del paso de Armijo. Es decir, cuando $N$ es mas grande entonces $h_k$ es mas pequeño. Ademas, consideramos un maximo de $100$ iteraciones.
%%time
problem_1.solve(
x0=[0, 2],
eps=1e-6,
N=10**2,
sgm=0.2,
verbose=False,
max_iter=100
)
----------------- Params ----------------- Start point : [0, 2] epsilon : 1e-06 (Armijo) N : 100 (Armijo) sigma : 0.2 (Armijo) max. iter.: 100 ------------------------------------------ --------------------------------------------------- Solver : Zoutendijk (simplex/Armijo) Solution time : 0.530 ms Objective : 48.0000 Stop criterion : max. iter. Max. iterations --------------------------------------------------- Wall time: 530 ms
Del reporte final obtenemos que el metodo tuvo que detenerse por exceder el maximo de iteraciones.
%%time
P1_report = problem_1.report(feasible=True)
P1_report
Wall time: 1 ms
| x1 | x2 | f(xk) | Factible | |
|---|---|---|---|---|
| Iter. k | ||||
| 0 | 0.000000 | 2.0 | 288.000000 | True |
| 1 | 2.000000 | 3.0 | 129.000000 | True |
| 2 | 4.000000 | 0.0 | 48.000000 | True |
| 3 | 3.333333 | 1.0 | 57.888889 | True |
| 4 | 4.000000 | 0.0 | 48.000000 | True |
| ... | ... | ... | ... | ... |
| 96 | 4.000000 | 0.0 | 48.000000 | True |
| 97 | 3.333333 | 1.0 | 57.888889 | True |
| 98 | 4.000000 | 0.0 | 48.000000 | True |
| 99 | 3.333333 | 1.0 | 57.888889 | True |
| 100 | 4.000000 | 0.0 | 48.000000 | True |
101 rows × 4 columns
En la siguiente parte calculamos el error entre el punto encontrado y la solucion, y el error entre los valores de la funcion objetivo en la ultima iteracion y el valor del problema.
%%time
problem_1.error(sol_P1, val_P1)
--------------------------------------------------- Last point : [4. 0.] Solution : [3.8583198 0.2125203] Objective : 48.0 Value : 46.902912 --------------------------------------------------- Wall time: 0 ns
{'error_sol': 0.25541761439707455, 'error_val': 1.0970876585634315}
Podemos observar que los errores son altos, verificando que el punto encontrado no es la solucion real
Dado que los resultados obtenidos no fueron los mejores, realizaremos un analisis de los hiperparametros del metodo para entender cual es una configuracion adecuada para resolver este problema. Para esto, solamente probaremos distintos valores de N y sigma. Para el primero se estudiaran valores de la forma $10^k$, con $k=1,...,5$. Y para el segundo, $21$ valores en $[1e6, 0.9]$. Con esto obtenemos una grilla de 105 pares de valores diferentes. Esto se traduce en 105 instancias diferentes de hiperparametros para resolver el problema 1.
%%time
grid = {
"N": [10**k for k in range(1, 6)],
"sgm": np.linspace(1e-6,0.9, 21)
}
print("Tamaño de la grilla:", len(grid["N"]) * len(grid["sgm"]))
problem_1_grid = Zoutendijk(
fun="8 * ( x1 - 6 ) ** 2 + ( x2 - 2 ) ** 4",
const=[
"- x1 + 2 * x2 <= 4",
"3 * x1 + 2 * x2 <= 12",
"x1 >= 0",
"x2 >= 0"
],
info=False
)
results_grid = []
for N in grid["N"]:
for sgm in grid["sgm"]:
o = {"N": N, "sgm": sgm}
problem_1_grid.solve(
x0=[0, 2],
eps=1e-6,
N=N,
sgm=sgm,
verbose=False,
max_iter=100
)
o = {**o, **problem_1_grid.error(sol_P1, val_P1)}
o = {**o, **problem_1_grid.status()}
o = {**o, **problem_1_grid.opt()}
results_grid.append(o)
Tamaño de la grilla: 105 Wall time: 28.8 s
df_grid = pd.DataFrame(results_grid)
df_grid[df_grid["feasible"]].sort_values("error_val", ascending=True)
| N | sgm | error_sol | error_val | dt | stop_criterion | n_iter | feasible | xf | f(xf) | |
|---|---|---|---|---|---|---|---|---|---|---|
| 95 | 100000 | 0.495000 | 2.709507e-07 | 1.158185e-12 | 0.121109 | lmbdk < eps | 11 | True | [3.8583196488781555, 0.21252052668276664] | 46.902912 |
| 96 | 100000 | 0.540000 | 7.722711e-07 | 9.400480e-12 | 0.103095 | lmbdk < eps | 9 | True | [3.858319370795645, 0.2125209438065317] | 46.902912 |
| 97 | 100000 | 0.585000 | 7.742804e-07 | 9.478640e-12 | 0.138124 | lmbdk < eps | 12 | True | [3.8583193696810962, 0.21252094547835373] | 46.902912 |
| 98 | 100000 | 0.630000 | 1.122385e-06 | 1.985256e-11 | 0.158144 | lmbdk < eps | 14 | True | [3.8583191765874885, 0.21252123511876583] | 46.902912 |
| 93 | 100000 | 0.405001 | 1.715075e-06 | 4.632028e-11 | 0.321292 | lmbdk < eps | 29 | True | [3.8583207505267585, 0.21251887420986054] | 46.902912 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 67 | 10000 | 0.180001 | 2.554176e-01 | 1.097088e+00 | 0.548499 | max. iter. | 100 | True | [3.9999999999999996, 0.0] | 48.000000 |
| 26 | 100 | 0.225001 | 2.554176e-01 | 1.097088e+00 | 0.507462 | max. iter. | 100 | True | [3.9999999999999996, 0.0] | 48.000000 |
| 28 | 100 | 0.315001 | 2.554176e-01 | 1.097088e+00 | 0.503420 | max. iter. | 100 | True | [3.9999999999999996, 0.0] | 48.000000 |
| 69 | 10000 | 0.270001 | 2.554176e-01 | 1.097088e+00 | 0.556509 | max. iter. | 100 | True | [3.9999999999999996, 0.0] | 48.000000 |
| 88 | 100000 | 0.180001 | 2.554176e-01 | 1.097088e+00 | 1.114013 | max. iter. | 100 | True | [3.9999999999999996, 0.0] | 48.000000 |
89 rows × 10 columns
Graficaremos la grilla y pintaremos las instancias segun el criterio de parada que utilizo el metodo para detenerse.
fig = px.scatter(df_grid, x="N", y="sgm", log_x=True, color="stop_criterion",
hover_data=df_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Criterios de parada para los valores de sigma y N"
)
fig.show(renderer="notebook")
Podemos observar que para valores de sigma menores que 0.4 la mayoria de las instancias se detienen por superar el maximo numero de iteraciones. Por otro lado, existe un porcentaje no menor de instancias que se detienen por obtener iteraciones infactibles. Para estos casos, no existe un patron claro. Por otro lado, el resto de instancias se detiene por un paso pequeño.
Graficamos solo las instancias que son factibles. Y se colorea la grilla segun el error obtenido comparando la funcion objetivo en la ultima iteracion con el valor real del problema 1.
subdf_grid = df_grid[df_grid["feasible"]]
fig = px.scatter(subdf_grid, x="N", y="sgm", log_x=True, color="error_val",
color_continuous_scale="jet_r", hover_data=subdf_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Error absoluto del valor del problema para los valores de sigma y N. <br>No se consideran aquellos casos que son infactible."
)
fig.show(renderer="notebook")
Se observa que en las mismas instancias en que el metodo se detuvo por superar el maximo de iteraciones, entonces tambien tiene un error alto.
De lo observado anteriormente, solo visualizaremos aquellas instancias que se detuvieron ya que el paso era muy pequeño. Esta vez, el eje x es el valor del problema, el eje y es el valor de sigma y se colorea cada instancia segun el N.
subdf_grid = df_grid[(df_grid["stop_criterion"] == "lmbdk < eps")]
subdf_grid["N"] = subdf_grid["N"].apply(str)
fig = px.scatter(subdf_grid, x="error_val", y="sgm", log_x=True, color="N",
hover_data=subdf_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Error absoluto del valor del problema (eje x) para los valores de sigma y N (en colores). <br>Solo consideran aquellos casos que se detienen por lmbdk < eps"
)
fig.show(renderer="notebook")
C:\Users\felip\AppData\Local\Temp/ipykernel_19640/368999457.py:2: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
Del grafico anterior, para casi cualquier valor de sigma, se observa que si se elige un N grande entonces, el metodo logra obtener soluciones mas cercanas a la solucion real.
Mismo grafico anterior, pero esta vez los colores son en funcion al tiempo de ejecucion.
fig = px.scatter(subdf_grid, x="error_val", y="sgm", log_x=True, color="dt",
color_continuous_scale='jet_r', hover_data=subdf_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Error absoluto del valor del problema (eje x) para los valores de sigma y N, con tiempo de <br>ejecucion en colores. Solo consideran aquellos casos que se detienen por lmbdk < eps"
)
fig.show(renderer="notebook")
Del grafico anterior no es claro si los parametros afectan mas que otros para obtener mejores soluciones en menor tiempo.
De la busqueda de hiperparametro realizada anteriormente, se eligen valores de N=$10^5$ y sigma igual a $0.495$. Para lo cual, se obtienen los siguientes resultados:
%%time
problem_1 = Zoutendijk(
fun="8 * ( x1 - 6 ) ** 2 + ( x2 - 2 ) ** 4",
const=[
"- x1 + 2 * x2 <= 4",
"3 * x1 + 2 * x2 <= 12",
"x1 >= 0",
"x2 >= 0"
],
info=True
)
problem_1.solve(
x0=[0, 2],
eps=1e-6,
N=10**5,
sgm=0.495,
verbose=False,
max_iter=100
)
P1_report = problem_1.report(feasible=True)
P1_report
-------------- Model Size -------------- Variables : 2 Constraints : 4 Eq. consts. : 0 Ineq. consts.: 4 ------------------------------------------ ----------------- Params ----------------- Start point : [0, 2] epsilon : 1e-06 (Armijo) N : 100000 (Armijo) sigma : 0.495 (Armijo) max. iter.: 100 ------------------------------------------ --------------------------------------------------- Solver : Zoutendijk (simplex/Armijo) Solution time : 0.135 ms Objective : 46.9029 Stop criterion : lmbdk < eps Successful solution --------------------------------------------------- Wall time: 137 ms
| x1 | x2 | f(xk) | Factible | |
|---|---|---|---|---|
| Iter. k | ||||
| 0 | 0.000000 | 2.000000 | 288.000000 | True |
| 1 | 2.000000 | 3.000000 | 129.000000 | True |
| 2 | 4.000000 | 0.000000 | 48.000000 | True |
| 3 | 3.584656 | 0.623016 | 50.266235 | True |
| 4 | 3.878482 | 0.182278 | 46.923897 | True |
| 5 | 3.821375 | 0.267937 | 46.971498 | True |
| 6 | 3.859008 | 0.211488 | 46.902937 | True |
| 7 | 3.857075 | 0.214388 | 46.902992 | True |
| 8 | 3.858333 | 0.212501 | 46.902912 | True |
| 9 | 3.858314 | 0.212529 | 46.902912 | True |
| 10 | 3.858320 | 0.212521 | 46.902912 | True |
| 11 | 3.858320 | 0.212521 | 46.902912 | True |
En solo 11 iteraciones el metodo logra llegar al optimo.
Para comparar nuestros resultados con el optimo real del problema $(P_2)$ consideramos el motor https://www.wolframalpha.com. Para esto, utilizamos una version simplificada, dado que dicho problema se puede transformar en uno unidimensional. Esto es:
$$(P_2) \quad \begin{cases} \min &x_1^4+10x_1^3-x_1^2-8x_1+16\\ s.a. &x_1\geq 1/2\\ &x_1\leq 4\\ &x_2=x_1\\ &x_3=2x_1-1\\ &x_4=4-x_1 \end{cases}$$La instancia de este problema se puede encontrar aqui. Que entrega la siguiente solución para $x_1$
| x1 | f(x1) |
|---|---|
| 0.5311288741492748261833066 | 13.04675363940572740382986 |
_x1 = 0.5311288741492748261833066
sol_P2 = np.array([_x1, _x1, 2*_x1-1, 4-_x1])
val_P2 = 13.04675363940572740382986
sol_P2
array([0.53112887, 0.53112887, 0.06225775, 3.46887113])
%%time
problem_2 = Zoutendijk(
fun="x1 ** 4 - 2 * x2 ** 2 + 10 * x1 * ( x2 ** 2 ) + x4 ** 2",
const=[
"x1 + x2 - x3 = 1",
"x1 + x4 = 4",
"x1 - x2 = 0",
"x1 >= 0",
"x2 >= 0",
"x3 >= 0",
"x4 >= 0"
],
info=True
)
-------------- Model Size -------------- Variables : 4 Constraints : 7 Eq. consts. : 3 Ineq. consts.: 4 ------------------------------------------ Wall time: 0 ns
Este es un problema de 4 variables, con 7 restricciones. 3 de igualdad y 4 de desigualdad.
Inicializamos con un punto incial, epsilon pequeño, sigma igual a $0.2$ y un $N=100$. Donde $N$ es tal que $h_k = \frac{\bar{\lambda}_k}{N}$ es el $h$ del paso de Armijo. Es decir, cuando $N$ es mas grande entonces $h_k$ es mas pequeño. Ademas, consideramos un maximo de $100$ iteraciones.
%%time
problem_2.solve(
x0=[2, 2, 3, 2],
eps=1e-6,
N=10**2,
sgm=0.2,
verbose=False,
max_iter=100
)
----------------- Params ----------------- Start point : [2, 2, 3, 2] epsilon : 1e-06 (Armijo) N : 100 (Armijo) sigma : 0.2 (Armijo) max. iter.: 100 ------------------------------------------ --------------------------------------------------- Solver : Zoutendijk (simplex/Armijo) Solution time : 0.878 ms Objective : 14.7103 Stop criterion : max. iter. Max. iterations --------------------------------------------------- Wall time: 878 ms
Del reporte final obtenemos que el metodo tuvo que detenerse por exceder el maximo de iteraciones.
%%time
P2_report = problem_2.report(feasible=True)
P2_report
Wall time: 1 ms
| x1 | x2 | x3 | x4 | f(xk) | Factible | |
|---|---|---|---|---|---|---|
| Iter. k | ||||||
| 0 | 2.000000 | 2.000000 | 3.000000 | 2.000000 | 92.000000 | True |
| 1 | 0.500000 | 0.500000 | 0.000000 | 3.500000 | 13.062500 | True |
| 2 | 0.818182 | 0.818182 | 0.636364 | 3.181818 | 14.710334 | True |
| 3 | 0.500000 | 0.500000 | 0.000000 | 3.500000 | 13.062500 | True |
| 4 | 0.818182 | 0.818182 | 0.636364 | 3.181818 | 14.710334 | True |
| ... | ... | ... | ... | ... | ... | ... |
| 96 | 0.818182 | 0.818182 | 0.636364 | 3.181818 | 14.710334 | True |
| 97 | 0.500000 | 0.500000 | 0.000000 | 3.500000 | 13.062500 | True |
| 98 | 0.818182 | 0.818182 | 0.636364 | 3.181818 | 14.710334 | True |
| 99 | 0.500000 | 0.500000 | 0.000000 | 3.500000 | 13.062500 | True |
| 100 | 0.818182 | 0.818182 | 0.636364 | 3.181818 | 14.710334 | True |
101 rows × 6 columns
En la siguiente parte calculamos el error entre el punto encontrado y la solucion, y el error entre los valores de la funcion objetivo en la ultima iteracion y el valor del problema.
%%time
problem_2.error(sol_P2, val_P2)
--------------------------------------------------- Last point : [0.81818182 0.81818182 0.63636364 3.18181818] Solution : [0.53112887 0.53112887 0.06225775 3.46887113] Objective : 14.7103 Value : 13.046754 --------------------------------------------------- Wall time: 1 ms
{'error_sol': 0.7594707030190523, 'error_val': 1.663580354173945}
Podemos observar que los errores son altos, verificando que el punto encontrado no es la solucion real
Al igual que los resultados obtenidos preliminarmente para el problema 1, los resultados obtenidos para el problema 2 no fueron los mejores. Por esta razon, realizaremos un analisis de los hiperparametros del metodo para entender cual es una configuracion adecuada para resolver este problema. Para esto, solamente probaremos distintos valores de N y sigma. Para el primero se estudiaran valores de la forma $10^k$, con $k=1,...,5$. Y para el segundo, $21$ valores en $[1e6, 0.9]$. Con esto obtenemos una grilla de 105 pares de valores diferentes. Esto se traduce en 105 instancias diferentes de hiperparametros para resolver el problema 2.
%%time
grid = {
"N": [10**k for k in range(1, 6)],
"sgm": np.linspace(1e-6,0.9, 21)
}
print("Tamaño de la grilla:", len(grid["N"]) * len(grid["sgm"]))
problem_2_grid = Zoutendijk(
fun="x1 ** 4 - 2 * x2 ** 2 + 10 * x1 * ( x2 ** 2 ) + x4 ** 2",
const=[
"x1 + x2 - x3 = 1",
"x1 + x4 = 4",
"x1 - x2 = 0",
"x1 >= 0",
"x2 >= 0",
"x3 >= 0",
"x4 >= 0"
],
info=False
)
results_grid = []
for N in grid["N"]:
for sgm in grid["sgm"]:
o = {"N": N, "sgm": sgm}
problem_2_grid.solve(
x0=[2, 2, 3, 2],
eps=1e-6,
N=N,
sgm=sgm,
verbose=False,
max_iter=100
)
o = {**o, **problem_2_grid.error(sol_P2, val_P2)}
o = {**o, **problem_2_grid.status()}
o = {**o, **problem_2_grid.opt()}
results_grid.append(o)
Tamaño de la grilla: 105 Wall time: 52.9 s
df_grid = pd.DataFrame(results_grid)
df_grid[df_grid["feasible"]].sort_values("error_val", ascending=True)
| N | sgm | error_sol | error_val | dt | stop_criterion | n_iter | feasible | xf | f(xf) | |
|---|---|---|---|---|---|---|---|---|---|---|
| 74 | 10000 | 0.495000 | 9.104678e-08 | 2.131628e-14 | 0.060055 | |f(xk) * dk| < eps | 6 | True | [0.5311288397368259, 0.5311288397368259, 0.062... | 13.046754 |
| 96 | 100000 | 0.540000 | 3.431488e-07 | 2.788880e-13 | 0.134122 | lmbdk < eps | 9 | True | [0.531129003847329, 0.531129003847329, 0.06225... | 13.046754 |
| 99 | 100000 | 0.675000 | 6.480101e-07 | 1.003642e-12 | 0.221203 | lmbdk < eps | 15 | True | [0.5311291190740681, 0.5311291190740681, 0.062... | 13.046754 |
| 98 | 100000 | 0.630000 | 6.739846e-07 | 1.078249e-12 | 0.194177 | lmbdk < eps | 13 | True | [0.5311291288915007, 0.5311291288915007, 0.062... | 13.046754 |
| 97 | 100000 | 0.585000 | 8.586687e-07 | 1.753264e-12 | 0.166152 | lmbdk < eps | 11 | True | [0.5311291986955287, 0.5311291986955287, 0.062... | 13.046754 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 85 | 100000 | 0.045001 | 1.008680e+00 | 3.109605e+00 | 1.483769 | max. iter. | 100 | True | [0.9123741237412374, 0.9123741237412374, 0.824... | 16.156359 |
| 21 | 100 | 0.000001 | 1.040081e+00 | 3.329883e+00 | 0.842775 | max. iter. | 100 | True | [0.9242424242424243, 0.9242424242424243, 0.848... | 16.376636 |
| 63 | 10000 | 0.000001 | 1.057677e+00 | 3.457235e+00 | 0.877798 | max. iter. | 100 | True | [0.9308930893089309, 0.9308930893089309, 0.861... | 16.503988 |
| 42 | 1000 | 0.000001 | 1.057777e+00 | 3.457967e+00 | 0.884805 | max. iter. | 100 | True | [0.9309309309309309, 0.9309309309309309, 0.861... | 16.504721 |
| 84 | 100000 | 0.000001 | 1.057852e+00 | 3.458517e+00 | 1.480211 | max. iter. | 100 | True | [0.9309593095930959, 0.9309593095930959, 0.861... | 16.505270 |
105 rows × 10 columns
Graficaremos la grilla y pintaremos las instancias segun el criterio de parada que utilizo el metodo para detenerse.
fig = px.scatter(df_grid, x="N", y="sgm", log_x=True, color="stop_criterion",
hover_data=df_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Criterios de parada para los valores de sigma y N"
)
fig.show(renderer="notebook")
Podemos observar que para valores de sigma menores que 0.4 la mayoria de las instancias se detienen por superar el maximo numero de iteraciones. Por otro lado, existe solo una instancia que se detienen por el criterio del gradiente producto la direccion es pequeño. Por otro lado, el resto de instancias se detiene por un paso pequeño.
Graficamos todas las instancias que son factibles y se colorea la grilla segun el error obtenido comparando la funcion objetivo en la ultima iteracion con el valor real del problema 2.
subdf_grid = df_grid[df_grid["feasible"]]
fig = px.scatter(subdf_grid, x="N", y="sgm", log_x=True, color="error_val",
color_continuous_scale="jet_r", hover_data=subdf_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Error absoluto del valor del problema para los valores de sigma y N"
)
fig.show(renderer="notebook")
Se observa que en las mismas instancias en que el metodo se detuvo por superar el maximo de iteraciones, entonces tambien tiene un error alto.
De lo observado anteriormente, solo visualizaremos aquellas instancias que se detuvieron ya que el paso era muy pequeño. Esta vez, el eje x es el valor del problema, el eje y es el valor de sigma y se colorea cada instancia segun el N.
subdf_grid = df_grid[(df_grid["stop_criterion"] != "max. iter.")]
subdf_grid["N"] = subdf_grid["N"].apply(str)
fig = px.scatter(subdf_grid, x="error_val", y="sgm", log_x=True, color="N",
hover_data=subdf_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Error absoluto del valor del problema (eje x) para los valores de sigma y N (en colores).<br>Solo consideran aquellos casos que se detienen por lmbdk < eps"
)
fig.show(renderer="notebook")
C:\Users\felip\AppData\Local\Temp/ipykernel_19640/4254695722.py:2: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
Del grafico anterior, para casi cualquier valor de sigma, se observa que si se elige un N grande entonces, el metodo logra obtener soluciones mas cercanas a la solucion real.
Mismo grafico anterior, pero esta vez los colores son en funcion al tiempo de ejecucion.
fig = px.scatter(subdf_grid, x="error_val", y="sgm", log_x=True, color="dt",
color_continuous_scale='jet_r', hover_data=subdf_grid.columns)
fig.update_traces(marker_size=10)
fig.update_layout(
title="Error absoluto del valor del problema (eje x) para los valores de sigma y N, con tiempo de <br>ejecucion en colores. Solo consideran aquellos casos que se detienen por lmbdk < eps"
)
fig.show(renderer="notebook")
Del grafico anterior no es claro si los parametros afectan mas que otros para obtener mejores soluciones en menor tiempo.
De la busqueda de hiperparametro realizada anteriormente, se eligen valores de N=$10^4$ y sigma igual a $0.495$. Para lo cual, se obtienen los siguientes resultados:
%%time
problem_2 = Zoutendijk(
fun="x1 ** 4 - 2 * x2 ** 2 + 10 * x1 * ( x2 ** 2 ) + x4 ** 2",
const=[
"x1 + x2 - x3 = 1",
"x1 + x4 = 4",
"x1 - x2 = 0",
"x1 >= 0",
"x2 >= 0",
"x3 >= 0",
"x4 >= 0"
],
info=True
)
problem_2.solve(
x0=[2, 2, 3, 2],
eps=1e-6,
N=10**4,
sgm=0.495,
verbose=False,
max_iter=100
)
P2_report = problem_2.report(feasible=True)
P2_report
-------------- Model Size -------------- Variables : 4 Constraints : 7 Eq. consts. : 3 Ineq. consts.: 4 ------------------------------------------ ----------------- Params ----------------- Start point : [2, 2, 3, 2] epsilon : 1e-06 (Armijo) N : 10000 (Armijo) sigma : 0.495 (Armijo) max. iter.: 100 ------------------------------------------ --------------------------------------------------- Solver : Zoutendijk (simplex/Armijo) Solution time : 0.067 ms Objective : 13.0468 Stop criterion : |f(xk) * dk| < eps Successful solution --------------------------------------------------- Wall time: 69.1 ms
| x1 | x2 | x3 | x4 | f(xk) | Factible | |
|---|---|---|---|---|---|---|
| Iter. k | ||||||
| 0 | 2.000000 | 2.000000 | 3.000000 | 2.000000 | 92.000000 | True |
| 1 | 0.500000 | 0.500000 | 0.000000 | 3.500000 | 13.062500 | True |
| 2 | 0.722272 | 0.722272 | 0.444544 | 3.277728 | 13.740221 | True |
| 3 | 0.541080 | 0.541080 | 0.082160 | 3.458920 | 13.048412 | True |
| 4 | 0.531068 | 0.531068 | 0.062136 | 3.468932 | 13.046754 | True |
| 5 | 0.531415 | 0.531415 | 0.062829 | 3.468585 | 13.046755 | True |
| 6 | 0.531129 | 0.531129 | 0.062258 | 3.468871 | 13.046754 | True |
En solo 6 iteraciones el metodo logra llegar al optimo.
En este trabajo nos hemos propuesto resolver en general problemas de optimización de tipo $(P)$ con restricciones lineales, tanto de desigualdad como de igualdad. Para resolver este tipo de problemas definimos el pseudoalgoritmo de las direcciones admisibles (o método de Zoutendijk).
Para implementar este método, las funciones que resuelven cada uno de los pasos del algoritmo se definen de forma independiente. Posteriormente, los pasos se ensamblan en una única función como método de una clase de Python, llamada $\texttt{Zoutendijk}$. Además de resolver problemas de tipo $(P)$, es capaz de leer automáticamente la función objetivo y las restricciones del problema a partir de un texto de estilo Python.
Se propusieron dos problemas de optimización para resolver con Zoutendijk. Cuyas soluciones se obtuvieron satisfactoriamente a partir de una búsqueda de hiperparámetros. Para cada problema se estudió el comportamiento del método para diferentes valores de N y sigma. Parámetros que son clave en el Paso de Armijo. El primero modula lo fino que es el dominio de búsqueda y el segundo caracteriza el dominio a partir de un recta. Para encontrar el mejor par de hiperparámetros para cada problema se utilizó la técnica GridSearch.
Del análisis de los parámetros obtuvimos algunas conclusiones interesantes. Entre estos estan: (1) Valores pequeños de sigma (menos de 0,4) hacen que el método utilice demasiadas iteraciones para alcanzar incluso soluciones no optimales, (2) Para valores mayores de sigma (mayores de 0. 4) el método consigue utilizar menos de 100 iteraciones y detenerse antes de tiempo debido a un paso pequeño, (3) Utilizar valores grandes de N (igual a 10^5 para el problema 1 y =10^4 para el problema 2) permite obtener soluciones más cercanas a la solución real, y (4) No se obtuvieron patrones claros para el tiempo de ejecución del código, sin embargo resolver cada problema no toma más de 1 segundo.
Finalmente, todo este trabajo está disponible en un repositorio de github. Véase el repositorio en github.com/furrutiav/nonlinear_optimization_zoutendijk. Además, quedan experimentos futuros para mejorar el rendimiento del método. Esto consiste en mejorar el método de Armijo con inicializaciones aleatorias de las variables y muestreo aleatorio del dominio de búsqueda. Por otro lado, queda por revisar el rendimiento del método para otros tipos de funciones objetivo que no sean polinómicas y regiones de decisión más intrincadas. Entre ellas se encuentra la conocida función Rosenbrock.